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JEFFREY HUMPHERYS, GREGORY LYNG, AND KEVIN ZUMBRUN 

Abstract. Extending recent results in the isentropic case, we use a 
combination of asymptotic ODE estimates and numerical Evans- function 
computations to examine the spectral stability of shock-wave solutions 
of the compressible Navier-Stokes equations with ideal gas equation of 
state. Our main results are that, in appropriately rescaled coordinates, 
the Evans function associated with the linearized operator about the 
wave (i) converges in the large-amplitude limit to the Evans function 
for a limiting shock profile of the same equations, for which internal 
energy vanishes at one endstate; and (ii) has no unstable (positive real 
part) zeros outside a uniform ball |A| < A. Thus, the rescaled eigenvalue 
ODE for the set of all shock waves, augmented with the (nonphysical) 
limiting case, form a compact family of boundary-value problems that 
can be conveniently investigated numerically. An extensive numerical 
Evans- function study yields one-dimensional spectral stability, indepen- 
dent of amplitude, for gas constant 7 in [1.2,3] and ratio v/fi of heat 
conduction to viscosity coefficient within [0.2,5] (7 « 1.4, u/fi ~ 1.47 
for air). Other values may be treated similarly but were not considered. 
The method of analysis extends also to the multi-dimensional case, a 
direction that we shall pursue in a future work. 
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1. Introduction 

A long-standing question in gas dynamics is the stability of viscous shock 
layers, or traveling- wave solutions 

U{x,t) = U(x- st), lim U(z) = U±, 

z— >±oo 

of the compressible Navier-Stokes equations, where U(x, t) = (y, u, E) T is 
a vector recording specific volume, velocity, and total energy of the fluid at 
location i£M and time t £ R + . A closely related question is the relation 
between Navier-Stokes solutions and solutions of the formally limiting Euler 
equations in the limit as viscosity and heat conduction coefficients go to 
zero: more precisely, validity of formal matched asymptotics predicting that 
the Navier-Stokes solution consists approximately of an Euler solution with 
smooth viscous shock layers replacing discontinuous Euler shocks. 

Recent progress in the form of "Lyapunov-type" theorems established in 
[381 154"! |22| [23| |24"] has reduced both problems to determination of spectral 
stability of shock layers, i.e., the study of the eigenvalue ODE associated 
with the linearized operator about the wave: a standard analytically and 
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numerically well-posed (boundary value) problem in ODE that can be at- 
tacked by the large body of techniques developed for asymptotic, exact, and 
numerical study of ODE. Indeed, the cited results hold for a substantially 
more general class of equations, and in one- or multi-dimensions. 

In [291 [44] [TBI I16j . it has been established in a similarly general con- 
text (general equations, one- and multi-dimensions), using asymptotic ODE 
techniques, that spectral stability holds always in the small-amplitude limit, 
where amplitude is measured by \U + — C7_|, i.e., for shocks sufficiently close 
to a constant solution, thus satisfactorily resolving the long-time stability 
and small- viscosity limit problems for small- variation solutions. 

However, until very recently, the spectral stability of large- amplitude 
shock waves has remained from a theoretical viewpoint essentially open, 
the sole exceptions being (i) a result of stability of Navier-Stokes shocks for 
isentropic gas dynamics with 7-law gas in the special case 7 — > 1, obtained by 
Matsumura-Nishihara [39J quite early on using an ingenious energy estimate 
specific to that case; and (ii) and a result of Zumbrun [55] — again obtained 
by energy estimates special to the model — which establishes the stability of 
stationary phase-transitional shocks of an isentropic viscous-capillary van 
der Waals model introduced by Slemrod [49J. 

Progress instead has focused, quite successfully, on the development of ef- 
ficient and general numerical methods for the study of stability of individual 
waves, or compact families of waves, of essentially arbitrary systems; see, 
for example, [6j [U El [30, 4J. These techniques, based on Evans-function 
computations, effectively resolve the question of spectral stability for waves 
of large but finite amplitude, but leave open the question of stability in the 
large- amplitude limit. For discussion of the Evans function and its numerical 
computation, see [H El EH El 11] or Section E3] below. 

Quite recently, however, Humpherys, Lafitte, and Zumbrun [27J have in- 
troduced a new strategy combining asymptotic ODE techniques with nu- 
merical Evans-function computations, by which they were able to carry out 
a global analysis of shock stability in the isentropic 7-law case, yielding 
stability independent of amplitude for 7 G [1,2. 5] a Specifically, after an 
appropriate rescaling, they showed by a detailed asymptotic analysis of the 
linearized eigenvalue ODE that the associated Evans functions (determining 
stability of the viscous shock profile) converge in the large-amplitude limit 
to an Evans function associated with their formal limit, which may then be 
studied either numerically or analytically (for example, by energy estimate 
as in [27]). 

The purpose of the present paper is to extend the approach of [27J to 
the full (nonisentropic) Navier-Stokes equations of compressible gas dynam- 
ics with ideal gas equation of state, resolving in this fundamental case the 
long-standing questions of viscous shock stability and behavior in the small- 
viscosity limit. Specifically, we show, as in the isentropic case, that the 
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Evans function indeed converges in the large-amplitude limit, to a value 
corresponding to the Evans function of a limiting system. Compactifying 
the parameter range by adjoining this limiting system, we then carry out 
systematic numerical Evans- function computations as in [U [27] to determine 
stability for gas constant 7 £ [1.2, 3] and (rescaled) ratio of heat conduction 
to viscosity coefficient v/fj, € [0.2, 5], well-including the physical values given 
in Appendices |A]and [HI The result, as in the isentropic case, is unconditional 
stability, independent of amplitude for an ideal gas equation of state. 

1.1. Discussion and open problems. The asymptotic analysis of [27] is 
quite delicate; it depends sensitively both on the use of Lagrangian coordi- 
nates and on the precise way of writing the eigenvalue ODE as a first-order 
system. It is thus not immediately clear that the analysis can be extended 
to the more complicated nonisentropic case. Moreover, since Lagrangian 
coordinates — specifically, the associated change of spatial variable 



where p is density — are not available in multi-dimensions, it is likewise, at 
first glance, unclear how how to extend the analysis beyond one spatial 
dimension. 

Remarkably, we find that the structure of the full, physical equations 
is much more favorable to the analysis than that of the isentropic model. 
In particular, whereas in the isentropic case the eigenvalue equations in the 
large-amplitude limit are a nonstandard singular perturbation of the limiting 
equations that must be analyzed "by hand" , in the full (nonisentropic) gas 
case, they are a regular perturbation for which convergence may be concluded 
by standard theorems on continuous dependence of the Evans function with 
respect to parameters; see, for example, the basic convergence lemma of [44] . 

Indeed, for 7 bounded away from the nonphysical case 7 = 1 (see Section 
[2] for a description of the equations and the physical background) , we have 
the striking difference that, for a fixed left endstate £/_, the density remains 
uniformly bounded above and below for all viscous shock profiles connecting 
U- to a right state U+, with energy going to infinity in the large-amplitude 
limit. By contrast, in the isentropic case, the density is artificially tied to 
energy and thus density goes to infinity in the large-amplitude limit for 
any 7 > 1; see, e.g., |46[ |4"7] 148] . This untangling of the large-amplitude 
behaviors of the density and the energy sets the stage for our analysis. 
Below, to see this untangling, instead of fixing a left endstate U- and asking 
which right endstates U+ may be connected to U- by a viscous shock profile, 
we fix the shock speed s = — 1 and all coordinates of the left state U- except 
the energy. We find again that the density stays bounded above and below 
for all possible right states J7+ connected by a shock profile to some such 



Since the equations remain regular so long as density is bounded from 
zero and infinity, one important consequence of this fact is that we need only 





U-. 
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check a few basic properties such as uniform decay of profiles and continuous 
extension of stable/unstable subspaces to conclude that the strong-shock 
limit is in the nonisentropic case a regular perturbation of the limiting system 
as claimed; see Section [3] for details. 

A second important consequence is that Lagrangian and Eulerian coor- 
dinates are essentially equivalent in the nonisentropic case so long as 7 re- 
mains uniformly bounded from 1, whereas, in the isentropic case, the equa- 
tions become singular for Eulerian (x) coordinates in the large-amplitude 
limit, by (jl.ip together with the fact that density goes to infinity. Here, 
we have chosen to work with Lagrangian coordinates for comparison with 
previous one-dimensional analyses in the isentropic case [U [2?] [12] . How- 
ever, we could just as well have worked in Eulerian coordinates, including 
full multi-dimensional effects, to obtain by the same arguments that the 
large-amplitude limit is a regular perturbation of the (unintegrated) limit- 
ing eigenvalue equation, and therefore the Evans functions converge in the 
limit also in this multi-dimensional setting. 

Likewise, uniform bounds on unstable eigenvalues may be obtained in 
multi-dimensions by adapting the asymptotic analysis of [22] similarly as 
we have adapted in Section [4] the asymptotic analysis of [37] . Thus, for 
7 uniformly bounded from 1, the analysis of this paper extends with suit- 
able modification to the multi- dimensional case, making possible the res- 
olution of multi-dimensional viscous stability by a systematic numerical 
Evans- function study as in the one-dimensional case. We shall carry out 
the multi-dimensional analysis in a following work [28]. 

Presumably, the same procedure of compactifying the parameter space 
after rescaling to bounded domain would work for any gas law with ap- 
propriate asymptotic behavior as p,e — * 00. Thus, we could in principle 
investigate also van der Waals gas/fluids, for example, which could yield 
interesting different behavior: in particular, (as known already from stabil- 
ity index considerations |57[ 154]) instability in some regimes. Other inter- 
esting areas for investigation include the study of boundary layer stability 
(see |12] for an analysis of the isentropic case), and stability of weak and 
strong detonation solutions (analogous to shock waves) of the compressible 
Navier-Stokes equations for a reacting gas. A further interesting direction 
is to investigate the effects of temperature dependence of viscosity and heat 
conduction on behavior for large amplitudes; see Appendices IB. 21 and El 

In this work, we have restricted to the parameter range 7 G [1-2,3] and 
y//i£ [0.2,5], where 7 is the gas constant, v — k/c v is a rescaled coefficient 
of heat conduction {k the Fourier conduction and c v specific heat), and (i 
is the coefficient of viscosity; see equations. (|2.ip ~ (|2.3p . Section [2] Similar 
computations may be carried out for arbitrary 7 bounded away from the 
nonphysical limit 7 = 1. To approach the singular limit 7 = 1 would 
presumably require a nonstandard singular perturbation analysis like that 
of [27] in the isentropic case, as the structure is similar; see Remark l2.2i The 
limits v/fjj — > and z///x ^ 00 are more standard singular perturbations with 
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fast / slow structure that should be treatable by the methods of PQ ; this would 
be a very interesting direction for further study. We note that our results 
for large u/fj, do indicate possible further simplification in behavior, as the 
singular perturbation structure would suggest; see Remark 14.81 and Figure 
HI For dry air at normal temperatures, 7 ~ 1.4 and u/fi ~ 1.47, well within 
range; see Appendix lAl 

Finally, we mention the issue of rigorous verification. Our results, though 
based on rigorous analysis, do not constitute numerical proof, and are not 
intended to. In particular, we do not use interval arithmetic. Nonetheless, 
the numerical evidence for stability appears overwhelming, particularly in 
view of the fact that the family {D(X, of Evans contours estimated in 
the stability computations is analytic in both parameters, yielding extremely 
strong interpolation estimates by the rigidity of analytic functions. 

In any case, our analysis contains all of the elements necessary for numer- 
ical proof, the effective realization of which, however, is a separate problem 
of independent interest. Given the fundamental nature of the problem, we 
view this as an important area for further investigation. 

1.2. Plan of the paper. In Section [21 we set up the problem, describing 
the equations, rescaling appropriately, and verifying existence and uniform 
decay of profiles independent of shock strength. In Section [3l we construct 
the Evans function and establish the key fact that it is continuous in all pa- 
rameters up to the strong-shock limit. In Sectional we carry out the main 
technical work of the paper, establishing an upper bound on the modulus 
of unstable eigenvalues of the linearized operator about the wave in terms 
of numerically approximable quantities associated with the traveling-wave 
profile. In Section we describe our numerical method, first estimating a 
maximal radius within which unstable eigenvalues are confined, then com- 
puting the winding number of the Evans function around the semicircle with 
that radius to estimate the number of unstable eigenvalues, for (a discretiza- 
tion of) all parameters within the compact parameter domain, including the 
strong-shock limit. Finally, in Section [6l we perform the numerical compu- 
tations indicating stability. 

In Appendices [A] and El we discuss further the dimensionless constants 
r and v/n, and determine their values for air and other common gases. In 
Appendix[Cl we discuss equations of state for fluids and dense gases. In Ap- 
pendix |Dj we compute a formula for the Mach number, a useful dimension- 
less quantity measuring shock strength independent of scaling. In Appendix 
IB"! we give a general bound on the operator norm of lifted matrices act- 
ing on exterior products, useful for analysis of the exterior product method 
of [H [5] . In Appendix El we discuss the changes needed to accommodate 
temperature-dependence in the coefficients of viscosity and heat conduction, 
as predicted by the kinetic theory of gases. 



SPECTRAL STABILITY OF IDEAL-GAS SHOCK LAYERS 



7 



2. Preliminaries 

In Lagrangian coordinates, the Navier-Stokes equations for compressible 
gas dynamics take the form 

(2.1) v t -u x = 0, 

(2.2) Ut+Px= (^ 

V v 

(2.3) + + (^k) , 

V V Jx \ v J x 

where v is the specific volume, u is the velocity, p is the pressure, and the 
energy E is made up of the internal energy e and the kinetic energy: 

u 2 

(2.4) E = e+ Y- 

The constants fi and k represent viscosity and heat conductivity. Finally, 
T is the temperature, and we assume that the internal energy e and the 
pressure p are known functions of the specific volume and the temperature: 

p = p (v,T), e = e (v,T). 

An important special case occurs when we consider an ideal, polytropic 
gas. In this case the energy and pressure functions take the specific form 

RT 

(2.5) Po(v,T) = , e (v,T) = c v T, 

v 

where R > and c v > are constants that characterize the gas. Alterna- 
tively, the pressure may be written as 

Te 

2.6 p=— , 

v 

where r = 7 — 1 = — >0,7> 1 the adiabatic index. Equivalently, in terms 
of the entropy and specific volume, the pressure reads 

p(v,S) = ae^^v-f, 

where S is thermodynamical entropy, or p(v) = av 1 in the isentropic ap- 
proximation; see [4"6"1 141 [27], 

In the thermodynamical rarefied gas approximation, 7 > 1 is the average 
over constituent particles of 7 = (N + 2)/N, where N is the number of 
internal degrees of freedom of an individual particle, or, for molecules with 
"tree" (as opposed to ring, or other more complicated) structure, 

(2.7) 7 = *i±», 

where n is the number of constituent atoms [3]: 7 = 5/3 ~ 1.66 for 
monatomic, 7 = 7/5 = 1.4 for diatomic gas. For fluids or dense gases, 
7 is typically determined phenomenologically [26J. In general, 7 is usually 
taken within 1 < 7 < 3 in models of gas-dynamical flow, whether phe- 
nomenological or derived by statistical mechanics [46l El [48] . 
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2.1. Viscous shock profiles. A viscous shock profile of (|2.ip - (|2.3p is a 
traveling-wave solution, 

(2.8) v(x, t) = v(x — st), u(x , t) = u(x — st) , T(x,t) = T(x — st), 

moving with speed s and connecting constant states (v±,u±,T±). Such a 
solution is a stationary solution of the system of PDEs 

(2.9) v t - sv x -u x = 0, 

(2.10) u t - su x + p Q (v,T) x = (— ^ ) , 

\ V / x 

(2.11) 



[ e o(^r)+n 2 /2] t - S [eo( V ,r)+n 2 /2] a; + (po(^T)n) :l . = (^) + 



2.2. Rescaled equations. Under the rescaling 

(2.12) (x, t, v, u, T) -> (-e«c, es 2 t,v/e, -u/(es), T/{e 2 s 2 )), 
where e is chosen so that V- = 1, the system (j2.9j) -( |27TTT) becomes 

(2.13) u t + - u x = 0, 

(2.14) u t + u x +p x 



jJ,U x 



v / x 



(2.15) /w • /•:, • inn),, (i^) 



\ v / x 



where the pressure and internal energy in the (new) rescaled variables are 
given by 

(2.16) p(y, T) = e- l s~ 2 p {ev, e 2 s 2 T) 
and 

(2.17) e{v, T) = e- 2 s- 2 e (ev, e 2 s 2 T); 

in the ideal gas case, the pressure and internal energy laws remain unchanged 

(2.18) p(v,T) = —, e(v,T)=c v T, 

v 

with the same R, c v . Likewise, T remains unchanged in (|2.6|) . 

2.3. Rescaled profile equations. Viscous shock profiles of (|2.13p - (|2.15|) 

satisfy the system of ordinary differential equations 

(2.19) v' - u = 0, 
(2-20) u ' +p{Vi T)'=(f^-\ , 

(2.21) [e(v,T) + u 2 /2]' + {p{v,T)u)' = f^-\ + f^Y , 
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together with the boundary conditions 

(u(±oo),u(±oo),T(±oo)) = (v±,u±,T±). 

Evidently, we can integrate each of the differential equations from — oo to x, 
and using the boundary conditions (in particular V— = 1 and u_ =0), we 
find, after some elementary manipulations, the profile equations: 



(2.22) 

fiv = v 
(2.23) 
kT' = v 



(v — 1) + p(v, T) — p(v- , T_ 



(v - l) 2 



+ v(v — 1) p(v-,T_) — (y — 1) 



We note that in the case of an ideal gas, with V- = 1, these ODEs simplify 
somewhat, to 



(2.24) 
(2.25) 



e'= V - 
v 



~[u(t>-l)+r(e-ue_)], 
A* 

^ V ~ l)2 + (e - e_) + (i; - l)re- 



where v := k/c v and T is as in (12. 6p . 

Remark 2.1. Remarkably, the right-hand sides of the profile ODE are poly- 
nomial in (v, e), so smooth even for values on the boundaries v = or e = 
of the physical parameter range. This is in sharp contrast to the isentropic 
case [H [27], f or which the ODE become singular as v — > 0, except in the 
special case 7 = 1. 

2.4. Rankine-Hugoniot conditions. Substituting u + , e+ into the rescaled 
profile equations (j2.19fl (|2.2in and requiring that the right-hand side vanish 
yields the Rankine-Hugoniot conditions 



(2.26) 
(2.27) 

(2.28) 



-s[v\ 
-s[u] 
u 2 

e + T 



-[p], 

-[pit], 



where [/(£/)] := f(U+) — f(U-) denotes jump between U±. 

We specialize now to the ideal gas case. Under the scaling (|2.12|) . we have 
s = — 1, V- = 1, u- = 0. Fixing T max > T > T m i n > and letting v + vary 
in the range 1 > v+ > v*(T) := T/(T + 2), we use P^ - O^gD to solve for 



u+, e + and e_ 
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Our assumptions reduce ()2.26|) - (I2.28|) to 

(2.29) v+ - 1 = u+, 

(2.30) u + = -(p + - P .) = -r( ei 



+ 

2 

2.31 e+ - e_ + -± = - P+u+ = -T^±. 

2 v + 

Equation ()2.29|) immediately gives u+. Subtracting ^y- times (I2.30P from 
(|2.3ip . and rearranging, we obtain 

1 + ^(1 -v+) e.v+ {T + 2-TV+)) 



(2.32) e+ = e_ 



2« + 

r+2 

(2.33) t>+ > : 



l-25z(l-*+) ( r + 2 ) 



= pxn, from which we obtain the physicality condition 



r + 2' 

corresponding to positivity of the denominator, with — — ► +oo as v — ► u*. 
Finally, substituting into 1 = s 2 = — M and rearranging, we obtain 

(2.34) e= (r + 2)K-„.) 



2r(r + 1) 

and thus 

(o M v + (T + 2-Tv + ) 

(2 - 35) e+ = 2 r ( r + 1) ■ 

completing the description of the endstates. 

We see from this analysis that the strong-shock limit corresponds, for 
fixed r, to the limit v+ — > u#, with all other parameters functions of «+. In 
this limit, 

(2.36) w_ = 1, u_ = 0, e_ ~ (u+ - u*) -> 0, 
and 

(2.37) u + ~ (v + - 1) -> — e + ' ~ 1 



r + 2' 2(r + i) (r + 2) 2 

At this point, taking without loss of generality fi = 1, we have re- 
duced to a three-parameter family of problems on compact parameter range, 
parametrized by T max >T> T min > 0, 1 > v + > v*(T) > v*(T min ) > 0, 
and v max > > v m i n . 

Remark 2.2. As T — ► 0, we find from ()2.35p i/tai e+ blows up as (?+ — v*)/T, 
i.e., our rescaled coordinates remain bounded only if v + — v* < CT — » 0, 
C > constant, as well. (This is reflected in the limiting profile equation 
for r = 0, which admits only profiles from V- = 1 to v+ = 0; see (|2.25p . 
which, for T = 0, reduces to v' = v(v — 1).) Thus, our techniques apply for 
r —* onZy in i/ie (simultaneous) large- amplitude limit. 
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2.5. Existence and decay of profiles. Specializing to the ideal g as case, 
we next study existence and behavior of profiles. Existence and exponential 
decay of profiles has been established by Gilbarg [21] for all finite- amplitude 
shocks 1 > v + > ?/*. Thus, the question is whether these properties extend to 
the strong-shock limit, the main issue being to establish uniform exponential 
decay as x — > ±oo, independent of shock strength. 

Since the profile equations (12.24p - (|2.25j) are smooth (polynomial) in (v , e), 
the issue of uniform decay reduces essentially to uniform hyperbolicity of end- 
states (v, e)±, i.e., nonexistence of purely imaginary linearized growth/decay 
rates at ±oo. Linearizing (|2.24p - (|2,25p about an equilibrium state, we ob- 
tain 

(,3S) Q'.mQ, M:^-; JU) ( 2 1 "_-„ 1 + - r r e e ; \ 

Since M is 2 x 2, its eigenvalues are 

trM ± VtrM 2 — 4 det M 
m = , 

and so hyperbolicity is equivalent to det M/0 and det M < or trM ^ 0. 
Computing, we have 



(2.39) det M = (v/(J,v) ^(T + 2)v - {T + 1)(1 + Te_)^ 
so that det M ^ is equivalent (for T > 0, hence v > v * > 0) to 

(2.40) (r + 2)«-(r + i)(i + r e _) ^o. 

At f = f_ = 1, this reduces to e_ 7^ r(r+i) ' or ' nsm S (|2-34|) to 

, + ' r + 2 >^-^ >o, 

except in the characteristic case u+ = 1, while 

trM = - Ye- + (u/fi)' 1 ) > ~ 2{V + +l) + ~ V > " 

At v = v+, (flOOj) reduces, using ([234]) . to 

detM = ( , + /H(^ + (r + 2) ^ + " 1) )< ' 
except in the characteristic case v+ = 1. Thus, for t>+ bounded from zero, 
hyperbolicity fails at x = ±00 only in the characteristic case V- = v+ = 1. 

Next, let us recall the existence proof of [21], which proceeds by the 
observation that isoclines v' = and e' = obtained by setting the right- 
hand sides of (|2.24j) and (12.25D to zero bound a convex lens-shaped region 
whose vertices are the unique equilibria U±, that is invariant under the 
forward flow of (I2.24l) - (|2.25j) . and into which enters the unstable manifold of 
£/_; recall that det M < at — 00, hence there is a one-dimensional unstable 
manifold. It follows that the unstable manifold must approach attractor £7+ 
as x — > +00, determining the unique connecting orbit describing the profile. 
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By the above-demonstrated hyperbolicity, this argument extends also to 
the case v+ = t>* (e_ = 0), yielding at once existence and uniform bound- 
edness of profiles across the whole parameter range (recall that (v, e)± are 
uniformly bounded, by the Rankine-Hugoniot analysis of the previous sec- 
tion), in particular for the limiting profile equations at t>+ = v* of 

(2.41) v' = - [v(v - 1) + Te] , 

(,4 2) 

Collecting facts, we have the following key result. 

Lemma 2.3. For T bounded and bounded away from the nonphysical limit 
r = 0, /i, v bounded and bounded from zero, and v + bounded away from the 
characteristic limit V- = 1, profiles U = (v,u,e) T of the rescaled equations 
(|2.24p ~ (|2.25p exist for all 1 > v+ > v*, decaying exponentially to their 
endstates U± as x — > ±oo, uniformly in r,t> + ,//, v. 

Proof. Existence, boundedness, and exponential decay of individual profiles 
follow from the discussion above. Uniform bounds follow by smooth depen- 
dence on parameters together with compactness of the parameter range. □ 

3. Evans function formulation 

From now on, we specialize to the ideal gas case, setting without loss of 
generality /x = 1. 

3.1. Linearized integrated eigenvalue equations. Defining integrated 
variables 

/x rx px 

vdy, u := / udy, E := / Edy, 
-oo J— oo J— oo 

we note that the rescaled equations (|2. 13|) — (|2.15p can be written in terms 
of the integrated variables in the form 

v t + v x - u x = 0, 

, ~ r (.e^ y) u X x 

(o-\\ U t + U x -\ ; — = -z— , 

{6.1) Vx Vx 

~ -2 ~ ~2 

Tu x {Ex-^f) u x u xx , u ( E x - ^f)x 



E t + E x + l — = + 

V x 

The integrated viscous profile, 



v := / vdy, u := udy, E := Edy, 

J —oo J —oo J —oo 

is a stationary solution of (|3.ip . Then we write 

v = v + v, u = u + u, E = E + E, 
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and we linearize (J3JJ) about the integrated profile. By an abuse of notation, 
we denote the perturbation by v, u, and E. Note also that the integrated 
profile always appears under an x-derivative; this explains the appearance 
of "hats" in place of "tilde-hats" in the expression below. Finally, we use 

-2 

the relationship E = e — \ to simplify some of the expressions. We obtain 

the linearized integrated equations 

(3.2) 

v t + v x - u x = 0, 



T(E X — uu x ) Te 
u t + u x -\ -^v x 



u xx 
V 



Ux 

tkV x , 



Tu(E x — uu x ) Te Teu uu xx u x uu x 
E t + E x H h — u x ^-v x = — h — u x —v x 



+ 



V V 

v{E x - uu x ) x 



v v 

Defining e := E — uu, subtracting u times the second equation from the 
third, and rearranging, we obtain, finally, the linearized integrated eigen- 
value equations: 



Xv + v 



0, 



Xu + u + — e H —u + 



(3.3) 

Ae + e' + 



Te ii x 
v 2 v 2 



„ i / u xx 


u + 


Te 


U x 




V 




V 



U a 





~ve x 


u' + 


V 2 



// 

-e . 



3.2. Expression as a first-order system. Following [3], we may express 
(|3.3p concisely first-order system 



(3.4) 




W = A(x, X)W, 










( ° 


1 








\ 




Xu v 


v~ l i) v~ 1 vii x — u xx 


\g(u) 


9<P) ~ 


h(U) 


(3.5) A(x,X) = 








X 


1 















1 






V o 


r Xv + Tu x 


Xv 


f(U) 


-X ) 


(3.6) 


w 


= (e,e',u,v,v') T , i = 


d 

dx ' 






where, using u x = 


v„- and (|2.24() with u = l. 









(3.7) 



g (U) := iz-^re-^ + l)^) 



i/ _1 re 



i/"Te 



v 

u + 1 



(v(v-l)+T(e-ve-)) 
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(3.8) 



f(U) 



u. 



Te 



+ {, = 2v - 1 - r e _ 



(3.9) h(U) :-- 



ex_ 

v 



,-i 



+ (e-e_) + (t)-l)re. 



Remark 3.1. Remarkably, similarly as for the profile equations, the entries 
of A are polynomial in (v,u,e). Thus, both profile and linearized eigenvalue 
equations are perfectly well-behaved for any compact subset ofT > 0. 

3.3. Consistent splitting. Denote by A±(X) := lim x ^± OQ A(x, X) the lim- 
iting coefficient matrices at x = ±oo. (These limits exist by exponential 
convergence of profiles U, Lemma [2T3l ) Denote by S± and U± the stable 
and unstable subspaces of A±. 

Definition 3.2. Following [1], we say that (|3.4j) exhibits consistent splitting 
on a given X- domain if A± are hyperbolic, with dim S+ and dimCL constant 
and summing to the dimension of the full space (in this case 5). 

By analytic dependence of A on A and standard matrix perturbation 
theory, S + and U- are analytic on any domain for which consistent splitting 
holds. 



Lemma 3.3. For all T > 0, 1 > v+ > v*, (|3.4p - (l3.5p exhibit consistent 
splitting on {5RA > 0} \ {0} ; with dim S+ = 3 and dim£7_ = 2. Moreover, 
subspaces S+ and U-, along with their associated spectral projections, extend 
analytically in X and continuously in T, u, V—, to {$tX > 0} for T > 0, v > 0, 
and 1 > v + > u*. 



Proof. Consistent splitting and analytic extension to A : 
general results of [37], except at the strong-shock limit 



(3.10) 



follow by the 
= v* , where 



and 



(3.11) 



A*+ (A) 







( 


1 





\ 
























(A) = 













A 1 



















1 










V 


r 


A 


A 1 - A/ 






/ 







i 











\ 


Xv 


-l 











5(^+ 


) 















X 


1 




















1 




\ 







r 


Aw* 


Aw* /([/+) 


-v 



where g(U+) = ii ^ 1 (w* - [y - 1)) and /([/+ 

The matrix A*_ is lower block-triangular, with diagonal blocks 

A 

A 



r-2 
r+2- 





Xv' 1 



1 

,-i 
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corresponding respectively to the scalar convection-diffusion equation 

@t "T c x ve xx 

and the isentropic case treated in [27] , The first has eigenvalues /i = 1=I=V ^ +4A , 
so satisfies consistent splitting on {5ftA > 0}\{0}, with analytic continuation 
(since eigenvalues remain separated) to 5ftA > 0. The second, as observed 

in [27], has eigenvalues fj, = — A, +4A , hence likewise satisfies consistent 
splitting on {5ft A > 0}\{0} and (since the single unstable eigenvalue remains 
separated from the two stable eigenvalues) continues analytically to 3?A > 0. 
Indeed, the unstable manifold has dimension two for all 5ftA > 0, hence is 
analytic on that domain. This verifies the proposition at x = — oo by direct 
calculation. 

At x = +oo, the computation is more difficult. Here, we refer instead to 
the abstract results of [37] . which assert that hyperbolic-parabolic systems of 
the type treated here, including the limiting case, at least for v* 7^ 0, exhibit 
consistent splitting on {5ftA > 0} \ {0}, with analytic extension to 5ftA > 0, 
so long as the shock is noncharacteristic, i.e., the flux Jacobian associated 
with the first-order part of the equations have nonvanishing determinants at 
x = ±00. These may be computed in any coordinates, in particular (v, u, e). 
Neglecting terms originating from diffusion, i.e., including only first-order 
terms from the left-hand side, we obtain from (|3.3p the flux Jacobian 

1 -1 

-Te/v 2 1 T/v 
Te/v 1 

which has determinant A = 1 — T 2 e/v 2 — Te/v 2 , giving for v+ = v* (e_ = 0) 
that A-oo = 1 > and, calculating at v + = v* that 

Te + /vl = 2, A +00 = -1 - 2T < 0. 

Thus, we may conclude by the general results of [37] that consistent splitting 
holds at both x = ±00 on {3f?A > 0} \ {0} for 1 > -u + > v*, with analytic 
extension to 5ft A > 0. □ 

Remark 3.4. We note that the results of [37] do not apply at x = —00, 
v + = v* , where e_ = leaves the physical domain. Specifically, at this value 
the genuine coupling condition of Kawashima |34] . hence the dissipativity 
condition of [37] fails, and so we cannot conclude consistent splitting; indeed, 
the eigenvalue fx = A (corresponding to the decoupled hyperbolic mode) is 
pure imaginary for any pure imaginary A. 

3.4. Construction of the Evans function. We now construct the Evans 
function associated with (|3.4p - (|3.5p . following the approach of [371 S3] • 

Lemma 3.5. There exist bases V~ = (Vf , V^)(\), V + = (V 3 + , F 4 + , F 5 + )(A) 
ofU-(X) and 5+ (A), extending analytically in A and continuously in T, u, v_ 
to {5ft A > 0} for T > 0, v > 0, and 1 > v + > v*, determined by Kato's ODE 

(3.12) V' = (PP - PP)V, 
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where P denotes the spectral projection onto S + , f/_, respectively, and ' 
denotes d/dX. 

Proof. This follows from Lemma 13.31 by a standard result of Kato [33J , valid 
on any simply connected set in A on which P remains analytic. □ 

Lemma 3.6. There exist bases 

W~ = {W{, W 2 -)(A), W + = (W+, W+, W+)(X) 

of the unstable manifold at — oo and the stable manifold at x = +oo of 
(|3.4j) - (13.51) . asymptotic to e A ~ x V~ and e A+x V + , respectively, as x — > =Foo ; 
and extending analytically in A and continuously in T, u, v+ to {5RA > 0} for 
r > 0, v > 0, and 1 > v + > v*. 

Proof. This follows, using the conjugation lemma of [40], by uniform expo- 
nential convergence of A to A± as x — > ±oo, Lemma 12.31 □ 

Definition 3.7. The Evans function associated with (|3.4p ~ (|3.5p is defined 
as 

(3.13) D(X) := det(W + ,W~)\ X=0 . 

Proposition 3.8. The Evans function D(-) is analytic in X and continuous 
in T,v,v+ on 5RA > and T > 0, u > 0, and 1 > v+ > v*. Moreover, on 
{3ftA > 0} \ {0} ; its zeros correspond in location and multiplicity with eigen- 
values of the integrated linearized operator C, or, equivalently with solutions 
of (I3.3P decaying at x = ±oo. 

Proof. The first statement follows by Lemma 13.61 the second by a stan- 
dard result of Gardner and Jones [El H9] , valid on the region of consistent 
splitting. □ 

Remark 3.9. Proposition \3.8\ includes in passing the key information that 
the Evans function converges in the strong-shock limit u + — ► v* to the Evans 
function for the limiting system at t> + = v*, uniformly on compact subsets 
o/{KA > 0}, as illustrated numerically in Fig. 

Remark 3.10. The specification in (|3.12p of initializing bases at infinity 
is optimal in that it minimizes "action" in a certain sense; see [31 j for 
further discussion. In particular, for any constant-coefficient system, the 
Evans function induced by Kato bases A3. 12H is identically constant in A. 
For, in this case, bases W + and W~ are given at x = by the values 
prescribed in (I3.12p . and both evolve according to the same ODE, hence W = 
(W~, W + ) x =o satisfies W = [P,P']W and D(X) := detW = constant by 
Abel's Theorem and the fact that tr[P, P'] = 0, where [P, P'] := PP' - P'P. 

Remark 3.11. More generally, det(V~ ,V + ) = constant in the "traveling 
pulse" case U+ = U-, by the argument of Remark \3.10\ whence the Evans 
function constructed here may be seen to agree with the "natural" Evans 
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function (independent of the choice ofV^) 

(3 - 14) E{X) - det(V+,V~) ~ V+-V- 

defined in |43j for that case. The latter may in turn be seen to agree with a 
(2-modified) characteristic Fredholm determinant of the associated linearized 
operator L about the wave |20| . formally equivalent to 

-Ut det 2 (L - A) 



det 2 (/+(L - X)- l (L-Lq)) 



det2(-Lo — A) ' 

where Lq denotes the (constant- coefficient) linearized operator about the 
background state U±. Our construction by Kato's ODE thus gives a nat- 
ural extension to the traveling-front case of the canonical constructions of 
[431 [20] in the traveling-pulse case, neither of which generalizes in obvious 
fashion to the traveling-front setting (the difficulty in both cases coming from 
the fact that det(V + , V-) may vanish). 

4. HlGH-FREQUENCY BOUNDS 

We now carry out the main technical work of the paper, establishing the 
following uniform bounds on the size of unstable eigenvalues. 



Proposition 4.1. Nonstable eigenvalues A of (|3.3p . i.e., eigenvalues with 
nonnegative real part, are confined for 7 > 1, < v+ < 1 to a finite region 
|A| < A, for any 



/ 1 T* I _i_ \ T* I _|_ 9 / T* J^* 

(4.1) A> 2max{l,i4max( 1 ~~ 1 1 ++ 'T /0 V ~ + — fx, A) 

x V v ' 

where 

(4.2) i^I(^ a ) : =E^xf^(-) 5 

^ are as defined in f|4.22|> below, and \ -\ is the matrix operator 
norm with respect to any specified norm on C 5 . 

Before establishing Proposition l4.ll we give a general discussion indicating 
the ideas behind the proof. For v + > v*, such high-frequency bounds have 
already been established by asymptotic ODE estimates in [37]. For v + = v*, 
the problem leaves the class studied in [37] (specifically, the dissipativity 
condition is neutrally violated as discussed in Remark I3.4p , hence requires 
further discussion. 

However, a brief examination reveals that the argument of [37J applies in 
this case almost unchanged. For, recall that the method of [37| to obtain 
high-frequency bounds was to decompose the flow of the first-order eigen- 
value equation for high frequencies into parabolic growth and decay modes of 
equal dimensions r = dim(u, e) with growth rates K/i ~ ilAj 1 / 2 , and hyper- 
bolic modes of dimension n — r, in the present case dimv = 1, with growth 



18 



HUMPHERYS, LYNG, AND ZUMBRUN 



rate ~ ±(3iA + 1) up to an exponentially decaying error term ~ e -6 ' 1 ', the 
delicate point being to separate decaying from growing hyperbolic modes. 

In the present, degenerate case, the hyperbolic rates are only ~ 3RA plus 
decaying error term, and so the final, delicate part of the argument in [37] 
does not apply. However, since there is only a single hyperbolic mode, 
this part of the argument is not needed. Specifically, the \\\ 1 / 2 /C spectral 
gap between parabolic and hyperbolic modes allow us for high frequencies to 
decompose the flow of the eigenvalue equation into the direct sum of growing 
parabolic modes blowing up exponentially at x = +oo, decaying parabolic 
modes blowing up exponentially at — oo, and a single hyperbolic mode that 
blows up exponentially at — oo for v + > v* and, though it does not blow up 
exponentially for t> + = v*, is in any case always transverse to the unstable 
manifold at x = — oo. 

To put things another way, the unstable manifold at x = — oo consists for 
|A| sufficiently large precisely of growing parabolic modes, which blow up 
at x = +oo. Thus, there exist no zeros of the Evans function, since these 
correspond to solutions belonging to both the unstable manifold at — oo and 
the stable (hence decaying) manifold at +oo. This shows the existence of 
uniform high-frequency bounds- it remains to establish quantitative bounds 
by keeping track of constants throughout the argument. 

Remark 4.2. A review of the above shows that the same argument applies 
whenever hyperbolic modes are uniformly decaying or growing, i.e., in the 
situation identified in [381 l54"t [22] that all hyperbolic characteristic speeds 
have a common sign. Likewise, the multidimensional case may be treated by 
essentially the same argument, following the generalization given in [22] . 



Proof of Proposition \4-l\ We carry out the argument in two steps. 

1. Preparation. Recasting (|3,4p . (|3.5p in the standard coordinates of 
[37] . as 



(4.3) 
(4.4) 

B(x,\) 



Z' = B(x,X)Z, Z = (v, u, e, u' , e') T 

( —A 





\(f(U)-v) Xv + Tu x 



1 
1 



f(U) 



o \ 


1 
r 

,-i„' 



\ Xh(U) v- x vu x -u xx \v~ l v g(U)-h{U) u^vj 

and we decompose B as B = \B\ + Bq with 

/ -1 \ 



(4.5) B 1 {x)= 

f(fi)-v v 

\ h(U) v- x v 0/ 
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/o 








1 


\ 











1 

















1 





Fu x 





f(U) 


r 


\o 


v~ x vu x - u xx 





90) ~ h(U) 





and 



(4.6) B (x) 



Noting that B\ is lower block-triangular, with (1 x 1) upper diagonal block 
— 1 strictly negative, and (4 x 4) lower diagonal block 

/0 0\ 

v 

\0 V~H 0y 

having all zero eigenvalues, we block-diagonalize by the lower block-diagonal 
transformation Z := TX, 



a 



1 

6 h 



( 



1 

-e h 
\ 



■= -(a + hY 






f(U)-v 
V HU) ) 



where, since a 2 



= 0, (J + a) 


- 1 =1 


— a, hence 








/ ^ 




( \ 


( 





\ 

















f(U)-v 


+ a 


f(U)-v 




-f(u) + 


V 


V HU) J 




V HU) J 


\ 


-h(U) 


J 



and 



rj-\ — lrpf 




& 



I 











o\ 
































—v x 














\m 














j{U) ■= y- 1 {{Te- -(v- l))v x + e x ^j , to obtain X' = CX, 



where 



(4.7) 



C = T~ X BT - T~ X T' = ACi + C , 



Ci(x,A) 



l-\ 





















































V 


















v~ l v 








/ 
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is in a variant of block Jordan form and 

(v - f(U) 

v-f(U) 
-h(U) 



(4.8) C (x,X) 



k(U) 



o 






o 1 o \ 

1 

1 

2/(17) -0 r 



v vu r — u rr 



g(u) 



where 
(4.9) 



k(U) := (2/(17) - v)(v - /(£/)) - Th(U) + v x 
lit}) := g(U)(v - f(U)) - u- l vh{U) - j(U). 
Making the further transformation X = Q~y, 



Q ■-- 



1 (3 
I 4 J ' 



Q- 



l -p 
o h 



(3~ A- 1 (0,0,l,0)(/ 4 + a)- 1 
= A- 1 (0,0,l,0)(/ 4 -a) 
= \-\-v, 0,1,0), 



(3' 






-v x 








o\ 















































v° 














we obtain y' = Dy, 



D = Q- l CQ - Q- l Q' 

= ADi + D + A _1 Z)_i + A- 2 /7_2, 



where 7>i = Ci, 
(4.10) 77 (a;,A) 
(4.11) 



(v-f(U) 
v-f(U) 
-h(U) 
k(U) 











1 



2f(U)-v 



o \ 


1 
r 



v l vu x - u xx g(U) 



V 1 v 



(m(U) 


-v 2 + f(U)v - Tu x + v x 





3(0 - f(U)) 


-r\ 





-v(v-f(U) 





v - f(U) 








-vh(U) 





-h(U) 








-vk(U) 





k(U) 





^ 


-vl(U) 





l(U) 
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(4.12) 



where 



D- 2 (x,X) 



/0 






—vm(U) 







m (U) :=v(v-f(U)) 

Making the "balancing" transformation y 

then obtain Z' = EZ, E = V^DV, where 
/-A 



m(U) 





k(U). 
= VZ, V 



o\ 







h 









(4.13) 



E = 






















X^v-H 





A l/2 







\ 




o / 



+ e, 



where 



(4.14) 



e = e + x~ 1/2 e„ 1/2 + A _1 e_i + A- 3/2 e_ 3/2 + A~ 2 e_ 



e (x) = 



(4.15) 9_ 1/2 ( 3 



^ - /(f>) 








\ 






« - /(f>) 















-h(U) 





















V{U) - 


v r 






I 





g(u) 


v~ x v J 






/ o 








3(v-f(U)) 


-r\ 











v-f(U) 














-h(U) 







k(U) 


Tu x 













\l(U) v- 


l vu x — 












(4.16) 0_!(x) = 



fm(U) 




V 



-v{v- f{U))+v x -Tu x 0\ 

-v(v-f(U) 

vh(U) 

k(U) 

l(U) oj 



(4.17) 



@-3/2(a 



/o 





Vo 



m{U) 0\ 





-vk(U) 

-vl(U) 0/ 
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(4.18) 



e_ 2 (x) 



(0 


—vm(U) 








0\ 















































\o 











0/ 



Finally, setting Z = SX, with 
'l 



s :- 



A := 



s) ' 
\fl 



J I 

-A A 



s- 1 



1 // 







2 V 7 ^ 
l/yfv 

o vv^^y ' 



S^S: 





S- 1 a, 



s s x 



we obtain 

(4.19) 

where 

(4.20) 
with 

(4.21) M + :-- 



(v x /4v) 
X' = (F + F)X, 
F = S^ES 

Al/2^1/2 



I —I 
-I I 



M_ 
M 4 





X^y/v) 1 / 2 )' 
M_ := 



-A 

-AV2£l/2 










-X^v/v) 1 / 2 , 



and 

(4.22) 

where 



F = s- 1 es-s- 1 s x 



(4.23) F (x) = 

(v-f(U) 
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(4.24) ^_ 1/2 (i 
/ 



-HU) 

-10) 
2^/i/u 

10} 
2^/lTjv 



-3y/$(v-f(U)) 
J " (v-f(U)) 

n(U) h(U)V{i 



2\f% 



2Vu~' 



^-%v-f(U)) 

n(U) , h(U)VH 



2^/v/v 



+ 



n(U) 



o 
o 






3Vi(v-f(U)) 

n(U) _ h(U)Vi 
2Vu~ 1 v 2 

g)f + # (v-f(U)) 

n0) _ h(U)VH 
2 



2yjv/i> 



(4.25) 



(m(U) 







q(U) 

-v(v-f0))+k0) 
2 

vh(U) l(U) 

2 + 2v / ^ rT 
-v(v-f0))-k0) 



vh(U) 
2 



10} 
2v^ 



q(U) 

-v(v-f0))-k(U) 
2 

vh0) _ 10) 
2 



2VV- 1 
-v(v~}0))+k(U) 
2 

vh0) 10) 
2 + 



2Vu- 








/ 



o\ 








(4.26) 



q(U) := -v(v -f{U)) -rOx+t), 



(4.27) 



(4.28) 



•^-3/2(3 







v° 

/o 







2. Tracking. Denoting X_ 








v l / 2 m{U) 


0\ 


ii l / 2 k{U) 
2 





v 1 / 2 k0) 















2v / ^ rT 


2v / ^ rT 


{) 1 / 2 /c(C/) 





-V 1 / 2 k0) 





2v / ^ rT 





v^ 2 l0) 
2V^ =T 




—vm(U) 





-vm(U) (A 











































0/ 





(Xi.Xa.Xa) 31 , X + = (X 4 ,X 5 f, and 



JF — ^4-4- 



we obtain from (j4TT9"j) - P~27J|) 

|X_|'< |^__||X_| + |.F_ + ||X h |, 



(4.29) 
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from which, defining £ := we obtain by a straightforward com- 

putation the Ricatti equation 

(4.30) C'< (-min{l,i/- 1 / 2 }{) 1 / 2 5RA 1 / 2 +|^__| + |^ ++ |)c+l^-+l + l^+-IC 2 - 
Denote by 

minil^- 1 ! 2 }^/ 2 ^ 1 / 2 - |.F__| - |^++| 

c±:= 

(4.31) r 



± 



{1, I/- 1 /2}^l/2S ftA l/2 _ |jp _| _ |^ ++ K2 \jr 



'mm 



the roots of 

(4.32) (-min{l,i/- 1 / 2 }{) 1 / 2 ^A 1/2 + |^_-_| + |^ ++ |)c+l^-+l + l^+-IC 2 = 0. 
Assuming for all x the condition 



(4.33) max{l, y i/2 } l^--l + l-^+l + VI^-+ll^+- [ < ^1/2 

C± are positive real and distinct, whence, consulting (|4.30p . we see that 
£' < on the interval < C < C+- 

It follows that f2_ := {( < £_} is an invariant region under the forward 
flow of (|4.19p ; moreover, this region is exponentially attracting for ( < ( + . 
A symmetric argument yields that f2+ := {£ > £+} is invariant under the 
backward flow of (|4.19|) . and exponentially attracting for ( > Spe- 
cializing these observations to the constant-coefficient limiting systems at 
x = — oo and x = +oo, we find that the invariant subspaces of the limiting 
coefficient matrices from which the Evans function is constructed must lie 
in f2_ and 0+, respectively. By forward (respectively, backward) invariance 
of f2_ (respectively, fi-f), under the full, variable-coefficient flow, we thus 
find that the manifold Span W~ of solutions initiated at x = — oo in the 
construction of the Evans function lies in f2_ for all x, while the manifold 
Span W + of solutions initiated at x = +oo lies in Q, + for all x. 

Since f2_ and f2+ are distinct, we may conclude that under condition 
(I4.33p . Span W + and Span W~ are transverse and the Evans function does 

not vanish. But gZQ implies (Q3l) . by KA 1/2 > together with QQgp . 

□ 

4.1. Universality and convergence in the high-frequency limit. The 

bounds obtaining from (|4.ip may in practice be rather conservative, as il- 
lustrated by the following example. 

Example 4.3. For the constant- coefficient scalar operator L := d 2 — ad x , 
write (L—X)w = as a first-order system W' = (A+@)W , W = (w, w' /\ l / 2 ) T , 

where A = A 1 / 2 ^ > ® = ^a) ' Block- diagonalizing A by W = 
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RZ, R = ^l)' we °^ a ^ n Z' = {A + @)Z, where A = diag(l, — 1), 

Q = R- l OR = (a/2) . Applying the analog of (fP3D on KA > 

0, where 5 = 25RA 1 / 2 > |A| 1//2 ; we obtain nonexistence of eigenvalues for 
2\a\, giving eigenbound |A| < 4|a| 2 . By contrast, standard elliptic 
energy estimates give |A| < |a| 2 /4 ; which, by direct Fourier transform com- 
putation, is optimal. Comparing, we see that the tracking bound is of the 
correct order, 0(\a\ 2 ), but with coefficient 16 times larger than optimal. 

This simple calculation may explain the ratio of roughly 10 between the 
nonisentropic bounds found by tracking in Section 15.21 and the isentropic 
energy bounds found in [U [27] . The following result may be used to gauge 
at a practical level the efficiency of the analytical tracking bounds by a 
(nonrigorous, but typically quite sharp) numerical convergence study. 



Proposition 4.4. On the nonstable half-plane 5RA > 0, 

lim 

|A|^oo 



(4.34) lim D(X)/e aXl/2 = constant, 



(4-35) 

«:=(! + v- 1 / 2 ) I f (v 1/2 (x) - v 1 ! 2 ) dx+ [ °°(v 1/2 (x) - v 1 / 2 ) dx ) real. 



o 

Proof. Reviewing the proof of Proposition 14.11 we And that the initial trans- 
formation T is asymptotically constant in A as |A| — > oo, and thus, the 
projection onto the "hyperbolic" mode corresponding to the 1-1 entry has 
the same property. It follows that the Kato ODE R' = (PP' - P'P)R used 
to propagate initializing bases at ±oo (see Section [5]), where P is the projec- 
tion onto the stable (respectively, unstable) subspace of A± , asymptotically 
decouples, yielding a constant (stable) hyperbolic basis element, and two 
stable and two unstable "parabolic" basis elements coming from the 4x4 
lower right-hand block of matrix E further below. But, the latter decouples 
into what may be recognized as a pair of first-order systems corresponding 
to the scalar variable-coefficient heat equations ut = u xx and ut = vu xx . Ex- 
plicit evaluation of the Kato ODE, similar to but simpler than the treatment 
of Burgers' equation in [27], Appendix D, then yields that the Evans func- 
tion for would be asymptotically constant if A(x, A) were identically equal 
to A^ for x < and A + for x > 0. See also Remark 13. 10[ which yields the 
same result in much greater generality. 

Though A is not constant for x ^ 0, the | A [-asymptotic flow may be de- 
veloped as in [37] in an asymptotic series in A -1 / 2 (respectively. A -1 ) in 
parabolic (respectively, hyperbolic) modes, to see that, up to an asymp- 
totically constant factor (coming from c(x) = 0(1) terms in eigenvalues of 

various modes, through e ^± oo '- c ^ -c '- = ' :00 ''-' d:c ), the flow from y to x is given 

asymptotically by e J « w and e J v y ' m parabolic 
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modes and e~ x ( y ~ x ^ in the hyperbolic mode (note: constant rate, so no re- 
sulting correction), whence, correcting for variation in the integrand from 
the constant-coefficient case, we obtain (j4.34f) . 

We omit the details, referring the reader to [27J 02] and especially [37] for 
similar but more difficult calculations. □ 

Remark 4.5. We see from (|4.34p that the asymptotic behavior of contours 
is independent of shock amplitude or model parameters, being determined up 
to rescaling of X by sgn a. This explains the "universal" quality of contour 
diagrams arising here and in \27\ [T2] . In practice, it is not necessary to com- 
pute a, since the knowledge that limit H4.34H exists allows us to determine 
a, C by curve fitting of log D(\) = logC + a\ l l 2 with respect to z := A 1 / 2 , 
for |A| S> 1. When D is initialized in the standard way on the real axis, so 
that D{\) = D{\), C is necessarily real. 

Remark 4.6. Restricting the limiting Evans function := Ce aX±/2 to the 
imaginary axis, A = it, t £ R, we obtain 

D\ir) = Ce a|r|1/2/2 (cos+isin)(±a|r| 1/2 ), 
predicting increasing winding about the origin as t — > oo. 

Since the limiting Evans function := Ce aAl/2 is nonvanishing, we may 
obtain practical high-frequency bounds by a convergence study onD^ D\ 
requiring, say, relative error < O.f to obtain a conservative but reasonably 
sharp radius. Here, may be estimated numerically using profile data 
in the formulae of Proposition 14.41 by curve-fitting as described in Remark 
14.51 or, more conventionally, by numerical extrapolation in the course of the 
convergence study. 

See, for example, the computation displayed in Figures[TH2l comparing the 
nonisentropic Evans function to its high-frequency limit Ce aX±/2 for T = 2/3 
and n = v = 1, at contour radii A = 25 and A = fO, respectively, where 
C and a have been determined by first taking limits along the positive 
real axis. Clearly, convergence in both cases has already occurred, whence 
radius A = 10 is sufficient to bound unstable eigenvalues, similarly as in the 
isentropic case [U [27] . By comparison (see Section [SJ , tracking estimates 
give the much more conservative bound A = 100.4. More extreme cases are 
depicted in Figures [3] and H] for /j, = 1, v = 5, and gas constants T = 2/3 
and r = 1/5, respectively. Note that convergence has already occurred at 
radius A = 40, which is thus sufficient to bound unstable eigenvalues; by 
contrast, the bounds obtained by tracking are A = 391.3 and A = 1755.6, 
respectively. 

These figures clearly indicate the universal behavior of the high-frequency 
limit. This can be seen also in the contour plots of Figure [6l comparing 
contours for the same model and radius at different values of v+. 

Remark 4.7. More generally, the argument of Proposition [J. 3J\ yields 
(4.36) L>(A)/e QAl/2+/3A = constant, 
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Figure 1. Universal behavior at high frequency: The im- 
ages of the semicircle of radius 25 under the Evans function 
and its universal approximant (|4.36j) (C, a, and [3 determined 
by curve fitting), for a monatomic gas, T = 2/3, /x = ^ = l, 
in the worst case t>+ = v* = 1/4. Agreement is nearly ex- 
act on the image of the outer, circular arc and most of the 
imaginary axis, with deviations for |A| small. 

a and (3 real constants, for all hyperbolic-parabolic parabolic system of the 
form studied in [371 EE] , where [5 corrects for variation in the rates of growth 
(respectively, decay) in hyperbolic modes, which are given to first order by 
A times their convection rates [37]. That (5 = in the present case is an 
accident of Lagrangian coordinates, for which hyperbolic modes are convected 
(in the rest frame of the shock) with the constant fluid velocity —s. In the 
more general case (|4.36p . asymptotic behavior of contours is determined (up 
to rescaling in X) by sgn a together with the additional parameter [3/a 2 . 

Remark 4.8. Figure [7] is particularly intriguing, showing convergence also 
for small frequencies. This suggests the conjecture that D* might converge 
identically to in the singular limit T — > 0, v j {i — > oo for all frequencies, 
small as well as large. This would be an interesting direction for further 
investigation. We remark that (i) the limit of D* as T — > is accessible by 
our techniques, Remark \2.S\ and (ii) the limit — * oo should be accessible 
by standard singular perturbation techniques. 

4.2. The small-amplitude limit. We mention in passing the following, 
related result noted in [27], regarding the small- amplitude limit v+ — ► 1. 
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Figure 2. Universal behavior: The images of the semicir- 
cle of radius 10 under the Evans function and its universal 
approximant (14.36[) . for T = 2/3, fj, = v = 1, u + = v* = 1/4. 
(Tracking radius = 100.4.) 

Proposition 4.9. The Evans function D converges uniformly as v+ — * 1 
on compact subsets of > 0} \ {0} ; r, u, fj, > to a constant C(T, v, /i). 

Proof. For |A| sufficiently large, this follows by Proposition 14. 41 together with 
the fact [H] that profiles in the small-amplitude converge to an approximate 
Burgers equation profile given by a symmetric tanh function, for which a 
vanishes to order |1— v+\ in (|4.35p . For |A| bounded, this follows as described 
in |27| by the fact that the Evans function also converges in the small- 
amplitude limit to the Evans function associated with the scalar, Burgers 
equation, which by direct calculation is constant. (The latter fact may 
be deduced, alternatively, by a simple scaling argument showing that, for 
Burgers equation, the small- amplitude limit and large-amplitude limits are 
equivalent.) □ 

The significance of Proposition 14.91 is that the exponential rate of decay of 
profiles to their endstates U± goes to zero in the characteristic limit v + — > 1, 
as seen in the proof of Lemma 12.31 Thus, we cannot immediately conclude 
as in the "regular" large-amplitude limit v + — > v* even that a limit exists 
as u+ — > 1. Moreover, a second consequence is that the computational do- 
main [— L_,L + ] on which we carry out numerical evaluation of the Evans 
function enlarges to \L±\ — * oo as u+ — > 1, since this must be taken roughly 
proportional to the inverse of the exponential rate for numerical accuracy; 
see Section 15.3.31 Thus, the boundary v + = 1 is not directly accessible by 
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Figure 3. Universal behavior: The images of the semicircle 
of radius 40 under the Evans function and its universal ap- 
proximant (j4.36|) . for T = 2/3, fi = 1, v = 5, = *y* = 1/4. 
(Tracking radius = 391.3.) 

numerical Evans function techniques, but requires a singular-perturbation 
analysis, either carried out numerically or else analytically as above. Alter- 
natively, small-amplitude instability may be ruled out by energy estimates 
as described in the introduction; see [39j [30l 13]. 

Remark 4.10. Similarly as described just below Remark \4-6\ for the high- 
frequency limit, the small- amplitude limit may be used to obtain a sharp but 
nonrigorou^ lower bound on the amplitude of unstable shocks by a conver- 
gence study requiring relative error < 1/2 between Evans values and their 
constant limit: or, still sharper, and more conveniently estimated, to re- 
quire relative error < 1/2 between the Evans function and its estimated 
high-frequency approximant. Convergence in the small- and large- amplitude 
limits it|_ — ► 1 and — > v* is illustrated numerically in Fig. 

5. Numerical protocol 

We now describe the numerical algorithm, based on approximate compu- 
tation of the Evans function, by which we shall locate any unstable eigen- 
values, if they exist, in our system, over the compact parameter range un- 
der investigation. Specifically, using analyticity of the Evans function in 
SRA > 0, we numerically compute the winding number around a large semi- 
circle B(0,A) n > 0} enclosing all possible unstable roots, obtaining 




Recall that rigorous small-amplitude bounds are available by energy estimates |30ll4]. 
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Figure 4. Universal behavior: The images of the semicircle 
of radius 40 under the Evans function and its universal ap- 
proximant (j4.36[) . for T = 1/5, fj, = 1, v = 5, v+ = = 1/11. 
(Tracking radius = 1755.6.). Note that the two curves are 
essentially indistinguishable, except that the universal ap- 
proximate does not loop but instead cusps near the origin. 



a count of the number of unstable eigenvalues within, and thus of the to- 
tal number of unstable eigenvalues. This approach was first used by Evans 
and Feroe [14] and has been advanced further in various directions in, for 
example, PI 021 El El El El ED] ■ 

5.1. Approximation of the profile. Following [U |27], we compute the 
traveling wave profile using MatLab's bvp4c routine, which is an adap- 
tive Lobatto quadrature scheme. These calculations are performed on a fi- 
nite computational domain [— L_,L + ] with projective boundary conditions 
M±(U — U±) = 0. The values of approximate plus and minus spatial infin- 
ity L± are determined experimentally by the requirement that the absolute 
error |J7(±L-|-) — U±\ be within a prescribed tolerance TOL = 10 -3 . This 
requirement is not too demanding in practice; we make more stringent re- 
quirements later in evaluating the Evans function. 

5.2. Bounds on unstable eigenvalues. We next estimate numerically the 
coefficients A) defined in (|4.2I) . (14. lj) . using the numerically generated 
profiles described above, generating an iterative sequence Aj+i := T(Aj), 

T(A) :=2max{l,i/}maxf l ^~ l + | /0 2v/J "-+ (x,A)Y. 

x \ v ' ' 
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385.3 
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1755.6 


0.4 


211.7 


182.3 


175.1 


325.0 


762.0 


0.6 


222.5 


123.4 


111.5 


198.4 


449.8 


2/3 


226.8 


114.3 


100.4 


175.3 


391.3 


0.8 


236.5 


103.9 


85.3 


142.6 


307.1 


1.0 


253.8 


100.8 


73.7 


113.8 


229.2 


1.2 


274.3 


106.6 


69.7 


98.7 


183.1 


1.4 


300.8 


117.5 


70.5 


91.6 


154.9 


1.6 


347.4 


131.8 


74.5 


89.8 


138.0 


1.8 


397.3 


148.7 


80.8 


91.8 


128.6 


2.0 


450.3 


167.5 


88.7 


96.7 


124.7 



Table 1. Tracking bound A* vs. T, v (i>+ = i>*). 



This is easily seen to converge, with odd terms monotone increasing and 
even terms monotone decreasing or the reverse, to a fixed point A* = T(A H< ), 
which, by Proposition 14.11 gives a bound |A| < A* on the maximum mod- 
ulus of unstable eigenvalues 3ft A > 0. Computations for a range of typical 
parameter values are displayed in Table [I] for the worst case v+ = t?*. Note 
the degradation of bounds for v S> // or v <C fj,, a consequence of multiple 
scales (stiffness). 

Remark 5.1. The poor rigorous tracking estimates obtained for v S> [i 
or v -C fj, could in principle be improved by a refined tracking estimate 
separating further the parabolic modes: that is, taking account of the presence 
of multiple parabolic scales; see [37], [44] for methodology. Ultimately, this 
should be treated by singular perturbation techniques as in pQ, separating 
out also fast/slow behavior of the background profile. 

5.2.1. High-frequency convergence study. Alternatively, we could obtain more 
realistic, but nonrigorous unstable eigenvalue bounds by a convergence study 
as |A| — ► oo, as described in Section [4.11 A convenient algorithm is to es- 
timate coefficients C, a of the high-frequency approximant D(\) ~ Ce Qv ^ 
by linear fit of log D ~ log C + a\/A as A goes to real positive infinity, 
then approximate by binary search the value A* at which relative error be- 
tween D(\) and Ce a ^ is less than TOL\ = W~ 1 on the positive semicircle 
|A| = A*, 5RA > 0, indicating convergence to this tolerance, hence nonva- 
nishing of D, for |A| > A*, 3ftA > 0. The resulting bounds are much less 
conservative than those obtained by rigorous tracking estimates of Section 
15.21 see for example the discussion following Remark 14.61 We do not pursue 
this here, as our main interest is in rigorous bounds. However, the ob- 
servation seems quite important for practical numerical testing, as typical 
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differences in radius are an order of magnitude or more, and the size of the 
radius appears to be the main limiting factor in computational efficiency. 

5.3. Approximation of the Evans function. We compute the Evans 
functions for comparison by two rather different techniques, both of which 
have been demonstrated to give good numerical results. 

5.3.1. The exterior product method. Following j6[ IU [5], we work with 
"lifted equations" 

W' = A {k) W, W := Wi A ••• A W fc , 

evolving subspaces encoded as exterior products of basis elements Wj , where 

(5.1) A (fc) (Wi A • • • A W k ) := (AW X A • • • A W k ) + ■ ■ ■ + (Wi A • • • A AW k ), 

defining W + and W~ as k + - and fc_ -products of bases {W- } and {W~} of 
the subspaces of solutions of W = AW decaying at +oo and — oo; in the 
present case, k + = 3, A;_ = 2. 

In this setting, the stable (respectively, unstable) manifold at +oo (re- 
spectively. — oo) corresponds to a single solution/vector, eliminating diffi- 
culties of "parasitic modes", etc.; see [HJ, [301 El S] for further discussion. We 
then compute the Evans function as 

D{\) = W + A W~| x= o 

or, alternatively, as D(\) = W + • W~| z =o, where W + is an appropriate 
solution of the adjoint equation W = -(A^)*W; see El El [301 H for 
further details. 

5.3.2. The polar coordinate method ("analytic orthogonalization" ) . An al- 
ternative method proposed in [30] is to encode W = jQ, where "angle" 
f2 = uj\ A • • • A ujk is the exterior product of an orthonormal basis {ujj} of 
Span {W±, . . . , Wfc} evolving independently of 7 by some implementation 
(e.g., Davey's method) of continuous orthogonalization and "radius" 7 is 
a complex scalar evolving by a scalar ODE slaved to 0, related to Abel's 
formula for evolution of a full Wronskian; see [30] for further details. This 
might be called "analytic orthogonalization", as the main difference from 
standard continuous orthogonalization routines is that it restores the im- 
portant property of analyticity of the Evans function by the introduction of 
the radial function 7 (f2 by itself is not analytic). 

Advantages/disadvantages: The exterior product method is linear, 
but evolving in a high-dimensional (~ n n ) space. The polar coordinate 
method is nonlinear, hence less well-conditioned and involving more compli- 
cated function calls, but evolving on a lower-dimensional manifold (~ n 2 ). 
Thus, there is a tradeoff in dimensions vs. conditioning, with the polar co- 
ordinate method the only reasonable option for high-dimensional systems 
and the exterior product method somewhat faster and more efficient for 
low-dimensional systems [30J. Both methods are effective (and reasonably 
comparable in efficiency) for n < 7 or so. 
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Role as numerical check: Since the two methods involve completely 
different ODE, one linear and the other nonlinear, agreement in their results 
is strong, if indirect, evidence that equations have been properly encoded 
and solutions accurately approximated, at least on the finite computational 
domain [— L_,L + ]. We address in the following subsection the separate 
question of determining appropriate L± . 



5.3.3. Determination of approximate spatial infinities. Denoting by A± (A) 
the limits at ±00 of the lifted matrix A^ k± \x,X) defined in (|5.ip . and //+ 
(respectively. \i- the eigenvalue of A^f} (respectively. A^) of smallest 
(respectively, largest) real part, we find that there holds a uniform bound 

(5.2) e (A«-M) ± * < a> 

on any compact subset of 5RA > 0, for T bounded from zero, and fx, v 
bounded and bounded from zero, for some C* > 0. For A bounded from 
zero, this follows by consistent splitting on {!RA > 0} \ {0}, and the choice of 
k± as dimensions of stable (respectively, unstable) subspaces of A±, which 

a (k) 

together imply that, away from A = 0, fi± are simple eigenvalues of A± . 
For A near zero, on the other hand, we may verify directly that ^± are 
semisimple, by the same considerations used to verify continuous extension 
of stable (respectively, unstable) subspaces of A±: simple in the case of fi-, 
since the unstable subspace remains uniformly spectrally separated from re- 
maining eigenvalues of A_\ semisimple in the case of //+, because hyperbolic 
characteristics a~j~ are simple, and thus eigenvalues fij of small real part, by 
the standard theory of |17[ l56l I37j . are analytic and semisimple, of form 
fij ~ —\/a~j, whence /i + (the sum of the k + = 3 eigenvalues of largest real 
part) is semisimple as well. 

Applying the "quantitative gap lemma" of Theorem C.2, [4], we have 
therefore that the relative error between the solution W ± (=bL-|-) at plus or 
minus approximate spatial infinity x = ±L± and the constant-coefficient 
initialization v= t e /x±= ' :i± is bounded by for 

(5.3) 6 := a\A {k H; A) - 4 fc) (A)lL 1{ ±L, ± oo)- 

Using the bound |ilf( fc )| < k\M\ established in Appendix [El and the asymp- 
totic behavior 

(5.4) Aj(x, A) - A jt ±{\) « Q je ~ 9±lxl , 



34 HUMPHERYS, LYNG, AND ZUMBRUN 

A =: XA\ + Aq for x — ► ±00, we may thus estimate 



k\A(±L ± ,X)-A ± (X)\ 



max|^ fe )(-, A) - ^i fc) (A)| L i (±Li±oo) < max k\A(; A) - A ± (X)\ L i (±L>±oo) 
|A|<A I A|<A 



~ max 

|A|<A 6± 

^ , \M±L±) - Aq,±\ + A|Ai(±£±) - Ai^j 
|A(±£±,0)-Afc(0)| 



\A(±L±,A) - At (A) - (A(±L±,0) - A±(0))| 



0± 

This gives a theoretical relative error bound of TOL (tolerance) between 
initializing basis at ±L and actual basis for the theoretical Evans function, 
provided that 

s \A(±L ± ,0) - A ± (0)\ + \A(±L ± ,A) - A±(A) - (A(±L±,0) - A±(0))K 

or, approximately, 
(5.5) 

|A(±L±,0) - A±(0)| + \A(±L±,A) - A±(A) - (A(±L ± ,0) - A±(0))| 

Remark 5.2. Alternatively, working directly from (|5.4|) . we may solve (|5.5p 
luzf/i equality for 

ft* t log + log A; + log(|Q± I + A|Qf |) + log + log TOL- 1 
1"-°) ^± ~ ^ • 

A reasonably good bound (noting insensitivity of log, and also cancellation 
in large A vs. small effects) for k = 2, A ~ 10 2 , TOL = 10~ 3 , |Qi| = 1 
for the sparse matrix Q±, C* = 10 2 , and throwing out 6 and \Qq\ terms as 
negligible for large |A|, is 

T Iog2 + 71ogl0 17 

— r ± — ™r ± 

Numerical algorithm Starting with the values needed for accurate pro- 
file approximation, we increment L± by some fixed step-size (here, we have 
chosen step-size 5) until (|53]) is satisfied, taking k = 2, TOL = 10" 3 , and 
conservatively estimating C* = 10 2 . In Table [2j we have displayed values of 
0±, L± computed from with TOL = 1(T 3 , C* = 10 2 , for various values 
of T and v + , with \x = 1 and v set to the value v = (3/4) 9 ^7 , 7 = T + 1, 
predicted by the kinetic theory of gases; see Appendix [Bj In principle, C* 
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could be estimated numerically for a more precise bound. In practice, con- 
vergence studies reveal these bounds to be rather conservative; see Table El 
or [H [27] in the isentropic case. 

Translation to x = 0. The above-described algorithm is designed to 
achieve relative accuracy TOL of W± at x = ±L+, whereas the accuracy of 
the Evans function is determined, rather, by their relative errors at x = 0. 
A complete description of the error must thus include also a discussion of 
possible magnification through the evolution from ±L+ to 0. However, as 
discussed in [U [7J, the flow toward x = is in fact quite stable, since, 
by construction, the modes W± we seek to approximate are the fastest 
decaying in the flow toward infinity, hence the fastest growing in backward 
flow toward zero, with all other modes decaying exponentially in relative 
size. Thus, in practice, the resulting magnification in error is negligible; see 
[3 [7] for further discussion or numerical studies. 

Remark 5.3. Since the polar coordinate method computes the same quantity 
W in different coordinates, the initialization error bounds derived for the 
exterior product method apply also for the polar coordinate method and so 
the same values L± may be used for both computations. 

5.3.4. ODE protocol. Following [3 El El 13 E3 H], to further improve the 
numerical efficiency and accuracy of the shooting scheme, we rescale W 
and W to remove exponential growth/decay at infinity, and thus eliminate 
potential problems with stiffness. Specifically, we let W(x) = x V(x), 
where fi~ is the growth rate of the unstable manifold at x = — do, and we 
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L 


rel(W + (0)) 


rel(W r -(0)) 


rel(D) 


20 


2.2(-2) 


5.8(-3) 


3.8(-3) 


25 


3.3(-3) 


8.8(-4) 


4.6(-4) 


30 


4.7(-4) 


1.3(-4) 


5.3(-5) 


35 


6.4(-5) 


1.8(-5) 


5.7(-6) 


40 


8.7(-6) 


2.6(-6) 


8.9(-7) 


45 


l.l(-6) 


3.5(-7) 


2.2(-7) 



Table 3. Convergence of W + (0), W~(0), and D as L in- 
creases from 25 to 50, incrementing by 5, for a diatomic gas 
(7 = 7/5 and v/fj, = 1.47. The values were computed for a 
quarter circle of radius R = 69 consisting of 90 points. Rel- 
ative errors were computed at each point and the maximum 
value along the contour is reported with the next higher value 
of L being the baseline. 



solve instead 

V'O) = (A( k \x,\)-ii-I)V(x). 

We initialize V(x) at x = —00 to be the eigenvector of A^(X) corresponding 

to . Similarly, it is straightforward to rescale and initialize W(x) at 
x = +00. To produce analytically varying Evans function output, the initial 
data V(— L_) and V(L+) must be chosen analytically using (13.121) . The 
algorithm of [8j works well for this purpose, as discussed further in [4l [30] . 

The ODE calculations for individual A are carried out using MatLab's 
ode45 routine, which is the adaptive 4th-order Runge-Kutta-Fehlberg method 
(RKF45), after scaling out the exponential decay rate of the constant- 
coefficient solution at spatial infinity, as described just above. This method 
is known to have excellent accuracy [5j [30] ; in addition, the adaptive refine- 
ment gives automatic error control. Typical runs involved roughly 300 mesh 
points per side, with error tolerance set to AbsTol = le-6 and RelTol = 
le-8. 

5.4. Winding number computation. We compute the winding number 
for fixed parameter values about the semicircle d{\ : 3ftA > 0, |A| < A} 
by varying values of A along 180 points of the contour, with mesh size 
taken quadratic in modulus to concentrate sample points near the origin 
where angles change more quickly, and summing the resulting changes in 
arg(D(A)), using QlogD(A) = argZ)(A)(mod27r), available in MatLab by 
direct function calls. As a check on winding number accuracy, we test a 
posteriori that the change in argument of D for each step is less than tt/8. 
Recall, by Rouche's Theorem, that accuracy is preserved so long as the 
argument varies by less than tt along each mesh interval. 
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Figure 5. Iso-Mach curves in V, v + . 



6. Numerical results 

Finally, we describe our numerical results in various cases, using the nu- 
merical method just described, varying t> + between .7 and the theoretical 
lower value v* in our rescaled coordinates (12.120 . For comparison between 
these and standard coordinates, it is useful to convert these values to their as- 
sociated Mach number, a standard dimensionless measure of shock strength 
discussed further in Appendix (Dj As computed in (|D. 1[) , this is given by 

M = ; = , which for 1 — v+ small reduces using ( 12.331 ) to 

1 2 + r 

M = = « 1 + (1 - v+)—— < 1 + 1 - v+ 

for the range T £ [0, 2] under consideration. In particular, for the upper 
limit v+ = .7 of our computations, we have on the reduced range < T < 1 
the exact upper bound M < (.55)" 1 / 2 , or approximately Mach 1.35, well 
within the small-amplitude regime. Recall that stability of small-amplitude 
shocks has been shown analytically in |29j . 

The smallest computed physical value v + — 1>* = 1CT 3 corresponds to Mach 
M ~ 100, at which we see already convergence of the Evans function to that 
of the nonphysical limit u + = «*, corresponding to Mach M = oo. For a 
visual comparison, we display iso-Mach (constant Mach number) curves in 
the T-v+ plane in Figure [5j 
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6.1. Description of experiments: broad range. In the main case con- 
sidered, of jj, and v of comparable but wide-ranging size, a first set of ex- 
periments was carried out in the range V S [-2,2], v £ [0.2,5], sampling at 
mesh points 

(r, v) e {0.2, 0.4., 0.6, 2/3, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0} 
x {0.2,0.5,1.0,2.0,5.0}, 

55 pairs in all, where, for each value (T,i/), v + was varied among 8 mesh 
points on a logarithmic scale from .7 to «*, for a total of 440 runs in all. In 
each of the cases that we examined, the winding number was zero, indicat- 
ing spectral stability and thereby nonlinear stability and existence of shock 
layers in the vanishing viscosity limit, by the results of [38} l54~l l22j l23"j 124] . 
These runs took 12 days to complete, of which 8 days were spent on the 
upper extreme case v = 5, and 2 days were spent on the lower extreme case 
v = .2, both out of physical range. 

6.2. Description of experiments: physical range. A second set of ex- 
periments was carried out for V values corresponding to monatomic, di- 
atomic, etc. gas on a refined mesh for v within the smaller, physical range 
v £ [1,2] indicated by Appendices lAl and lB| sampling at 

(I» e {2/11,2/9,2/7,2/5,2/3} 

x {1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0}, 

with f + again varied among 8 mesh points on a logarithmic scale from .7 to 
t>*, again a total of 440 runs. The results again were winding number zero 
for each case tested, indicating spectral stability. These runs took 10 days 
to complete. We remark that runs for v = 5 and v = .2 took over twice as 
long to complete compared with the average, again reflecting the stiffness 
alluded to in Remark 15. 11 associated with difference in scale between v and 
H = l. 

Appendix A. Gas constants for air 

In this appendix, we list the various relevant gas constants for dry air at 
normal temperatures and pressures. For the specific gas constant (universal 
gas constant R 8.3 mol e S . K divided by molar mass), we have 

287.05 — ] — , 
kg-K 

J = Joules, kg = Kilograms, K = degrees Kelvin [ID] . For specific heat at 
constant volume, we have 

c„ 716 



kg-K 

at 0°C (degrees Celsius) [50] . Computing the dimensionless quantity Y 
R/c v , we thus obtain 

r « 0.401, 
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1.5 



0.5- 



E 0- 



-0.5 




0.5 



1.5 



Re 



Figure 6. Convergence to the limiting Evans function as 
v+ — > v* for a monatomic gas, V = 2/3, \i = v = 1, v* = 1/4. 
The contours depicted, going from inner to outer, are images 
of the semicircle of radius 25 under the Evans function D for 
v+ = .7, .6, .5, .4, .35, .3, .27, .26, .255, .251, .2501. The outer- 
most contour is the image under the limiting Evans function 
D* , which is essentially indistinguishable from the images for 
v+ = .251 and v+ = .2501. 



in remarkable agreement with the value T = .4 predicted by statistical me- 
chanical approximation T = 2 n+i ^ or a diatomic gas, n = 2. (Recall that 
air is composed of roughly 78% nitrogen, 21% oxygen, and 1% neon, so it is 
essentially a diatomic mixture.) 

For thermal conductivity, the ratio of heat flux to temperature gradient 
asserted by Fourier's law of heat conduction, we have the value 

W 

K ~ .025 -, 

m • K 

W = Watts = Joules per second, and for dynamic or "first" viscosity, the 
rate of proportionality of shear stress to velocity gradient of a shear flow 
asserted in Newton's law of viscosity, the value 



Hi w (1.78) x 10 



-5 k g 



m • s 

m = meters, s = seconds [TOj, [52]. Collecting these values, we may compute 
the "const ant- volume Prandtl number" 

(716)(1.78) x 10" 5 



(A.l) 



Pr„ 



(.025) 



.510, 
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or very nearly 1/2. This is related to the usual (constant-pressure) Prandtl 
number Pr := Cpfii/n by Prt, = Pr/7, where 7 := c p /c v is the heat capacity 
ratio, or adiabatic index, relating specific heat at constant pressure c p to 
specific heat at constant volume c v , under ideal gas assumptions, 7 = 1 + T. 
For a variety of gases over a fairly wide range of temperatures [52], p. 43], 
Pr 0.7. 

Recall that the effective viscosity appearing in the one-dimensional Navier- 
Stokes equations is /J, = 2/ii + /X2, the sum of twice the dynamic viscosity and 
the "second viscosity" fi2, which is commonly taken as \ii = — (2/3)/xi based 
on the assumption that pressure equals "mean pressure" defined as one-third 
the trace of the three-dimensional stress tensor, an approximation that ap- 
pears to be in good agreement with experiment at least for monatomic and 
diatomic gases [HJHS]. We therefore take \x ~ (4/3)/ii, giving 

(A.2) v//jt = (3/4)?^ « 1.47. 

Interesting values for computation are thus T = .4 (7 = 1.4), v/fj, = 1.47, 
modeling air. 

Remark A.l. A convenient alternative formula involving commonly tabu- 
lated dimensionless quantities is 

(A.3) u/fi = (3/4) 7 Pr-\ 

where Pr denotes (usual) Prandtl number and 7 = c p /c v the adiabatic index. 
Assuming the value Pr ~ 0.7 typical of simple (e.g., monatomic, diatomic) 
gases, we obtain the general rule of thumb 

(A.4) u/n « (1.09)7; 

see Table\4\for more precise values. 

Remark A.2. We note also the alternative formula 



C-i, 



-1 



■2-1 
.R 

for the Gruneisen constant in terms of specific heat at constant pressure. 
(For air, Cp « U/g • K = 1, 000 J/kg • K.J 

Appendix B. Gas constants for other gases 

B.l. Ideal Gases. Using the dimensionless formula (|A.3p . we next investi- 
gate typical parameter values for some common gases. For an ideal gas, the 
Prandtl number is predicted by Eucken's formula [9] 

(B.l) p r = ^_ 

giving in combination with (1A.3|) the simple relation 

(B.2) ^ = (3 /4)^f^. 
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Gas 


7 (th.) 


Pr (th.) 


Pr (exp.) 


v/y. (th.) 


u/fi (exp.) 


rel. err. 


He 


5/3 


.667 


.694 


1.88 


1.80 


4% 


Ar 






.669 




1.87 


.05% 


H 2 


7/5 


.737 


.712 


1.43 


1.47 


2% 


N 2 






.735 




1.43 


0% 


2 






.732 




1.43 


0% 


CO 






.763 




1.38 


3.5% 


co 2 


4/3 


.762 


.819 


1.31 


1.22 


6.8% 


H 2 S 






.929 




1.08 


17.6% 


S0 2 






.833 




1.17 


10.7% 


CH 4 


5/4 


.8 


.777 


1.17 


1.28 


9.5% 



Table 4. Theoretical vs. experimental values of Pr and 
is/fx and relative error in v/fj, at 20 °C = 68 °F (room tem- 
perature). 



We display in Table 0] the relation for various gases between the theoretical 
value (|B.1|) and experimentally measured values for Pr (Table 1.9-1 of [9j, 
adapted from [35], p. 250). Though not displayed, experimental values 
of the adiabatic index 7 = c p /c v match theoretical predictions to within 
1% for monatomic and diatomic gases, 3% for triatomic, and 4.8% for five- 
atomic gas CH4, as do experimental values vs. theoretical predictions of the 
Gruneisen coefficient T = R/c v . 

In summary, the equation of state and temperature law predicted by ideal 
gas theory appear to be extremely accurate at usual temperatures, while the 
predictions involving viscosity (//, Pr, k) degrade with molecular complex- 
ity, being nearly exact for monatomic gases, quite good for most diatomic 
gases, but only a first approximation for triatomic and more complicated 
gases. Indeed, the derivation of viscosity and heat conduction formulae 
involves additional simplifying assumptions whose validity degrades with 
structural complexity [9]. Thus, we may use with some confidence the theo- 
retical prediction (IB, 2ft for simple gases, but must consult experimental data 
for complex gases. 

Remark B.l. From (|2.7p and (jB.ip we obtain the theoretical range 

7 g [1,1.66...], u/fie [.75,1.88]. 

B.2. Temperature dependence and the kinetic theory of gases. 

Though the ratio (jB.lj) of viscosity and heat conductivity predicted by the 
kinetic theory of gases is constant, the absolute size predicted for either 
one depends on temperature, T. For example, the predicted viscosity for a 
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monatomic gas obtained through Chapman-Enskog expansion of the Boltz- 
mann equations with hard-sphere potential is Chapman's formula 



(B.3) Ml = Ml ( T ) = ( 2 /3)W— - 

V 7T(T Z 

where m is mass per particle in kilograms kg, k is the Boltzmann constant 
in Joules per degrees Kelvin J/K, T is temperature in degrees Kelvin K, and 
a is the collision cross section in meters squared m 2 , given approximately 
by one-half diameter squared j3j [9]. This appears to give good physical 
agreement, and a refined version given by Sutherland's formula 

(l + m)T 3 / 2 

(B.4) mi = ' , 

1 + m 

rh constant, to give extraordinarily good agreement [9]. More generally, 
viscosity is typically modeled by 

(B.5) /ii = CT 9 , 1/2 < g < 1, 

with q ~ .76 for air [36J. 

Properly, we should include the above-described temperature-dependence 
in the physical study of large-amplitude shock layers. Though beyond the 
scope of the present project, this appears to be feasible by a slight modifi- 
cation of the same techniques, as we discuss further in Appendix [FJ 



Appendix C. Liquids and dense gases 

For comparison, values of the Prandtl number Pr for various media at 
20 °C are pj p. 80]: 

• around 0.024 for mercury, 

• around 0.7 for air and helium, 

• around 7 for water, 

• around 16 for ethyl alcohol, 

• around 10, 000 for castor oil, and 

• around 12, 000 for glycerin. 

The adiabatic index (specific heat ratio) of water is 7 ~ 1.01 ~ 1.0, hence 

u/fj, = 7/Pr « Pr" 1 w .143, 

quite far from the gas values of Table [H 

Moreover, for dense gases or liquids, the ideal gas assumptions break 
down, and the polytropic equation of state (|2.6[) must be replaced by more 
sophisticated versions such as Peng-Robinson or "stiffened polytropic" equa- 
tions of state [261 E5]- F° r example, water is often modeled by a stiffened 
equation of state 



(C.l) 



P = Tpe - 7P0 
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behaving like a prestressed material, with base stress Pq and T determined 
empirically: for example, 7 ~ 6.1 and 

P = 2, 000 MPa = 2 x 10 9 N/m 2 = 2 x 10 9 kg/m • s 2 

or 7 = 7.42 and Pq = 296.2 MPa |32l I25j . Similar techniques are used to 
model liquid argon, nickel, mercury, etc. [26j ITT]. 

It would be very interesting to investigate the effects on stability of these 
modifications in the polytropic equation of state. For the moment, what we 
can say, physically, is that insofar as they conform to a polytropic (ideal gas) 
equation of state, shock waves are stable. However, above a critical density, 
even standard, simple (e.g., monatomic or diatomic) gases are observed not 
to conform to a polytropic law [13] . and in this regime our mathematical 
conclusions hold no sway. 

Appendix D. Computation of the Mach number 

A dimensionless quantity measuring shock strength is the Mach number 

u_ — a 



M 

c_ 

(for a left-moving shock), where U- is the downwind velocity, a is the shock 
speed, and c_ is the downwind sound speed, all in Eulerian coordinates. 
The conservation of mass equation in Eulerian coordinates is pt + (pu) x = 0, 
giving jump condition o~[p] = [pu], or, in Lagrangian coordinates, 

U + V- — U-V + 

a = . 

V- — v + 

Thus, 

u- — a v-(u- — u + ) V-[u] V- 
C_ C-(V- — v + ) c-[v] c_ ' 

which, under our normalization s = —1, V- = 1, gives 



(T+2)(v+— u») 
or, using e_ = 2 r(r+i) 



M = cz l = (r(i + r) e _)- 1 / 2 

V2 



(D.l) M = ^ 

y/(T + 2)(v+-v*) 

In the strong-shock limit v + — ► v*, M ~ (v+ — f*) _1//2 ; in the weak shock 
limit v+ —* 1, M 1. 

Appendix E. Lifted Matrix bounds 

We establish the following useful bound on the operator norm of the 
"lifted" matrix A^ k ' acting on /c-exterior products V = V\ A • • • A Vf- by the 
operation 

A {k) V := AVi A • • • A V k + ■ ■ ■ + Vy A • • • A AV k . 



44 



HUMPHERYS, LYNG, AND ZUMBRUN 



induced by a given matrix A, where by convention is represented with 
respect to standard basis elements ej 1 A • • • ej k . 

Lemma E.l. In the matrix operator norm \ ■ \ p with respect to £ p , 
(E.l) \A {k) \ p < k\A\ p , l<p<oo. 

Proof. It is sufficient to establish (jE.ip for p = 1 and p = oo, the result for 
other p following from the Riesz-Thorin interpolation theorem 

\M\ p <\M\ e s \\M\% 

9j > 0, 6i + #2 = 1, for si < p < s 2 . 

The i l operator norm is equal to the maximum over all columns of the 
sum of moduli of column entries, or, equivalently, the maximum i 1 norm of 
the image of a standard basis element. Applying this definition to A^ k ), we 
find readily that the £ norm of the image of a standard basis element is 
bounded by the sum of the i norms of k terms of form 

Ae h A-Ae ife , 

expanded in standard exterior product basis elements. As each of these 
clearly has i l norm bounded by the i 1 norm of Aej 1 , and thus by \A\\, we 
obtain the result for p = 1. The result for p = oo may be obtained by 
duality, using (A^)* = (A*)^ together with |M|oo = |M*|i. □ 

Appendix F. Temperature-dependent viscosity 

Finally, we discuss the changes needed to accommodate a common tem- 
perature or other dependence in the coefficients of viscosity and heat con- 
duction. For concreteness, we focus on the case 

(F.l) n(e) = /i*e 9 , «(e) = is*e q , 0<q<l, 

constant, encompassing the Chapman formula (|F3.3[) predicted by the 
kinetic theory of gases as well as the more general formula (|B.5P , indicating 
afterward by a few brief remarks the extension to more general situations. 

F.l. Rescaling. Under (j2.5l) . (IF.lj) . it is easily checked that the Navier- 
Stokes equations are invariant under the modified rescaling 

(F.2) 0, t, v, u, T) (-es\es\~ 2g x, es 2 \es\~ 2g t,v/e, -u/(es),T/(e 2 s 2 )), 

consisting of (|2.12p with x and t rescaled by the common factor | es | 2,3 . For 
the Chapman viscosity q = 1/2, this reduces to 

(F.3) (x, t, v, u, T) — > (— (sgn s)x, \s\t,v/e, —uj (es), T / (e 2 s 2 )), 

essentially undoing the rescaling in the x-coordinate. Evidently, the Rankine- 
Hugoniot analysis of Section 12.41 goes through unchanged. 
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F.2. Profile equations. The profile equations (|2.24p ~ (|2,25p are unaffected 
by dependence of n, k. Setting := k*/c v , and making the change of 
independent variable 

dx 

(F.4) — = n = ^e q , 

dy 

we may thus reduce them to the form 
(F.5) v' = — [v(v - 1) + r(e - ve-)] , 

(v - l) 2 



(F.6) e' = - 



+ (e - e_) + (v - l)Te. 



already treated in the constant- viscosity case, ' = To obtain the profile 
in original x-coordinates, we have only to recover the change of coordinates 
x = x(y) by solving (|F.4[) with e = e(y), e the constant-viscosity profile. 

Recalling that e{y) = e± + 0{e~ e \ y \), 6 > for y ^ 0, where e+ > 
and e_ > except in the strong-shock limit v+ — * v*, we find that x and 
y are equivalent coordinates on x > 0, and on x < are equivalent for v + 
bounded from the strong-shock limit u*. However, in the strong-shock limit 
v + = v*, e_ = 0, we have the interesting phenomenon that the x < part 
of the shock profile terminates at a finite value X_ = x(—oo), since 



-oo p— oo 



\x(-oo)\ = I / (e(y)/c v ) q dy\ < C / e"^ 1 dy < +oo. 
Jo Jo 

Remark F.l. Holding (v,u,e)- fixed, and varying v + toward its minimum 
value , we find that 

as p + ~ e + ~ (e+/e_) ~ (u+ — since u + is bounded from zero, e_ 

is being held constant, and ratio e+ /e_ is independent of scaling so may be 
computed from formulae (I2.34p - (I2,35p . Thus, shock width in the constant- 
viscosity case is of order \v+ — v*\ going to zero as shock strength (measured 
in specific volume v ) goes to its maximum value of\l — | . By comparison, 
for the Chapman viscosity fx ~ e 1 / 2 , the shock width remains approximately 
constant in the strong-shock limit, a significant deviation in the theories. 
See, e.g., [5H [26] for further discussion of this and related points. 

Remark F.2. Clearly, the same procedure may be used to determine profiles 
for arbitrary jj, = fj,(v,u,e), k//j, constant, setting ^ = fj,(v,u,e) in (IF.4p . 

F.2.1. Limiting behavior. The limiting profile equations at v+ = u* are 
(F.7) v = — [v(v - 1) + re] , 

l) 2 

(F.8) e' = — - K - 1 +e 
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Linearizing about U- gives 
(F.9) ( V \ ~ M 



M :-- 







-i 



,-i 



For > 1, we have, therefore, that the slow unstable manifold at — oo 

is tangent to (r(i/*//i*)/(l — (^*//x*)), 1), with growth rate ~ e~ u * y , and 
thus generically 

1 Uy _Vy T/^* 

e e X — v^j^ 



(F.10) 



rN_- ■ 



as y 



-oo. 



For ^//i* < 1, the slow manifold is tangent to (1,0), and we have the 
opposite situation that, generically, ^ ~ ^f- — * ±oo as y — » — oo. 

F.3. Eigenvalue equations. Computing the linearized integrated eigen- 
value equations as in Section 13.11 and making the change of variables (|F,4I) , 
we obtain after some rearrangement, the modified, temperature-dependent 
equations 
(F.ll) 

Xfxv + v — u = 0, 



Xfxu + 



1 + 



A/k + e' + 



ge y _ qu y - 
ev jlev. 

V*Uyy 



U 



u + 



u + 



r _ quy 

v ev 



€ + 



-u + 



V 2 V 2 



V 



Te u y qu^u y ey 
— - (i/* + 1)^- 



+ 



it 



'* // 



' = where (v,u, e) = (w, u, e){y) are just as in the constant- viscosity case. 
This yields a first-order eigenvalue system 



(F.12) 
(F.13) 

A(y,X) 

(F.14) 
where 

(F.15) 



W' = A(y,X)W, 



/ 1 

\p,v~H v~ l v v~ l vu y 





y q y — t}Uy 



u 



mi 





Xp,g(U) 
Xfi 




\ 

g(U)-h(U) 
1 
1 



Xfiv + Tu y X(jlv f-) f(U) - A/t / 

„ i\T 



f(U) : = ,+ 



2v 







1 = 






_ gu| 


et) 


(lev 






e 


/ie / 
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(F.16) g(U) := (re - (y. + l)u y ) - 

(F.17) h{U) := -| = -v^ (-^^ + (e - e_) + (t> - l)r e _) , 
(F.18) A == e 9 , 

reducing for q = to that of the constant-viscosity case. We omit the details, 
which are straightforward but tedious. 

Noting that all terms not appearing in the constant-viscosity case involve 
derivatives of the profile, hence vanish at y = ±00 so long as e± (appearing 
in denominators) do not vanish, we may conclude by the constant-viscosity 
analysis consistent splitting away from the strong-shock limit v+ — > u*, 
e_ — > 0. We may thus construct an Evans function that is analytic in A and 
continuous in all parameters away from the strong-shock limit. 

F.3.1. Limiting behavior. Assume, as for the physical cases considered above, 
that 

(F.19) u/fj, = u«//jl* > 1. 

Then, by the limiting analysis (IF.lOh . all terms in A(y, A) remain bounded 
and well-defined in the strong-shock limit. Thus, we may hope for conver- 
gence as in the constant-viscosity case. 

On the other hand, new terms e y /e, u y /e of order 




exhibit singular behavior in the v+ — ► v*, e_ — > limit reminiscent of that 
of the isentropic case [27 . In particular, since e_/e approaches its limiting 
value 1 as y — > — 00 only as |e — e_|/e_, this means that \A — A±\ does not 
decay at uniform exponential rate as y — > —00, but only as O^lV^), 
6 > 0, so that the strong-shock limit is a singular perturbation in the sense 
of \44\ [27| and not a regular perturbation as in the constant-viscosity case. 
Thus, an analysis of behavior in the strong-shock limit is likely to involve a 
delicate boundary-layer analysis similar to that of [27] in the isentropic case. 
This appears to be a very interesting direction for further investigation. 

Remark F.3. The appearance of condition (|F.19|) is unexpected, dividing 
into two cases the physical behavior in the strong-shock limit. 

Remark F.4. Our numerical experiments, though still quite preliminary 
(restricted to diatomic gas, T = A, = 1.43, q = 0.5) so far again 

indicate unconditional stability, also in the temperature- dependent case. 
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F.4. General dependence. We conclude by discussing briefly the case of 
general, possibly inhomogeneous dependence of viscosity on (v,u,T). The 
homogeneous case goes exactly as before, working in y-coordinates and not- 
ing Remark IF. 21 and the discussion of Section IF. 31 

The inhomogeneous case is well-illustrated by Sutherland's formula (lB.4j) . 
Fixing a left-hand state (the "true" strong-shock limit), without loss of 
generality V- = 1, and taking v + — > rescale by (|F,2j) with q = 1/2 and 
e = 1, s = s(v+). The result in the rescaled coordinates is 

F.20 fi = 4/3K U >\ 1 = 4/3 1 } _ 2 _ , 

[s z l ) + m 1 + s z m 

converging in the strong-shock limit v + — > v*, s ~ |u+ — — ► oo to the 
homogeneous Chapman formula 

H = (4/3)(l +m)T 1 / 2 . 

Other inhomogeneous laws go similarly, converging in the strong-shock limit 
in each case to an appropriate homogeneous law. Thus, inhomogeneous 
dependence poses no essential new difficulty. 
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